##User defined parameters and functions

print(params)
## $readParams
## [1] TRUE
# set to true if want to run for a limited number of rows (i.e. for code testing)
readParams <- params$readParams

makePredictions <- function(predictionDF, modelObject) {
  ##
  # predictionDF <- climDatPred
  # modelObject <- bestLambdaMod_grassShrub_totalHerb
  # ##
  
  # pretend to scale the input variables so the model object can predict accurately
  predictionDF <- predictionDF %>% 
  mutate(across(all_of(prednames), base::scale,scale = FALSE, center = FALSE)) 
  
  # modelPredictions
  modelPreds <- predict(object = modelObject, newdata = predictionDF, type = "response")
  # add predictions back into the input data.frame
  predictionDF <- predictionDF %>% 
    cbind(modelPreds)
  
  # truncate all predictions to max out at 100 
  #predictionDF[predictionDF$modelPreds>100 & !is.na(predictionDF$modelPreds),"modelPreds"] <- 100
predictionDF[predictionDF$modelPreds < 0 & !is.na(predictionDF$modelPreds),"modelPreds"] <- 0

  # print predicted data
 return(predictionDF)
}

Dependencies

library(tidyverse)
library(sf)
library(terra)
library(kableExtra)
library(knitr)
library(USA.state.boundaries)
library(tidyterra)
library(ggpubr)
library(USA.state.boundaries)

Read in necessary data

# Load Data ---------------------------------------------------------------
# data ready for modeling (that has been scaled)
modDat_1_s <- readRDS("./models/scaledModelInputData.rds")
  
# get the soil raster, which we'll use for rasterizing the imagery
soilRastTemp <- readRDS("../../../Data_processed/SoilsRaster.rds") %>% 
terra::unwrap()

# make a map of the predictions
test_rast <-  rast("../../../Data_raw/dayMet/rawMonthlyData/orders/70e0da02b9d2d6e8faa8c97d211f3546/Daymet_Monthly_V4R1/data/daymet_v4_prcp_monttl_na_1980.tif") %>% 
  terra::aggregate(fact = 12, fun = "mean")

# download map info for visualization
data(state_boundaries_wgs84) 

cropped_states <- suppressMessages(state_boundaries_wgs84 %>%
  dplyr::filter(NAME!="Hawaii") %>%
  dplyr::filter(NAME!="Alaska") %>%
  dplyr::filter(NAME!="Puerto Rico") %>%
  dplyr::filter(NAME!="American Samoa") %>%
  dplyr::filter(NAME!="Guam") %>%
  dplyr::filter(NAME!="Commonwealth of the Northern Mariana Islands") %>%
  dplyr::filter(NAME!="United States Virgin Islands") %>%
  dplyr::filter(NAME!= "U.S. Virgin Islands") %>% 
  sf::st_sf() %>%
  sf::st_transform(sf::st_crs(test_rast))) 
  #sf::st_crop(sf::st_bbox(modDat_1_sf)+c(-1,-1,1,1))

## add ecoregion boundaries (for our ecoregion level model)
regions <- sf::st_read(dsn = "../../../Data_raw/Level2Ecoregions/", layer = "NA_CEC_Eco_Level2") 
## Reading layer `NA_CEC_Eco_Level2' from data source 
##   `/Users/astears/Documents/Dropbox_static/Work/NAU_USGS_postdoc/cleanPED/PED_vegClimModels/Data_raw/Level2Ecoregions' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2261 features and 8 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -4334052 ymin: -3313739 xmax: 3324076 ymax: 4267265
## Projected CRS: Sphere_ARC_INFO_Lambert_Azimuthal_Equal_Area
regions <- regions %>% 
  st_transform(crs = st_crs(test_rast)) %>% 
  st_make_valid() 
ecoregionLU <- data.frame("NA_L1NAME" = sort(unique(regions$NA_L1NAME)), 
                        "newRegion" = c(NA, "Forest", "dryShrubGrass", 
                                        "dryShrubGrass", "Forest", "dryShrubGrass",
                                       "dryShrubGrass", "Forest", "Forest", 
                                       "dryShrubGrass", "Forest", "Forest", 
                                       "Forest", "Forest", "dryShrubGrass", 
                                       NA
                                        ))
goodRegions <- regions %>% 
  left_join(ecoregionLU)
mapRegions <- goodRegions %>% 
  filter(!is.na(newRegion)) %>% 
  group_by(newRegion) %>% 
  summarise(geometry = sf::st_union(geometry)) %>% 
  ungroup() %>% 
  st_simplify(dTolerance = 1000)

## contemporary climate data
climDat_temp <- readRDS("/Users/astears/Documents/Dropbox_static/Work/NAU_USGS_postdoc/PED_vegClimModels/Data_processed/EcoRegion_climSoilData.rds")

climDat_timeTemp <- climDat_temp %>% 
  mutate(spatialID = paste0(round(x,4), "_", round(y,4))) %>% 
  drop_na()

spatID_lu <- data.frame("spatialID" = unique(climDat_timeTemp$spatialID), 
                        "uniqueID" = 1:length(unique(climDat_timeTemp$spatialID)))

climDat_timeTemp <- climDat_timeTemp %>% 
  left_join(spatID_lu)

## because of the spatial joining we did, there are some "locations" that have data for less than the true time span we want (2010-2020), or multiple observations per year... remove those so we have the even representation over the time span we're looking for 
set.seed("12011993")
goodTable <- data.frame(table(climDat_timeTemp$uniqueID))
goodIDs_temp <- as.numeric(goodTable[goodTable$Freq == 11,"Var1"])

#goodIDs <- unique(names(table(climDat_timeTemp$uniqueID))[(table(climDat_timeTemp$uniqueID)==11)])
# select 1000 random IDs
goodIDs <- goodIDs_temp[sample(x = 1:length(goodIDs_temp), size = 1000, replace = FALSE)]

## get climate data for 500 different points randomly chosen across CONUS for all years in the contemporary climate/weather/soils dataset 
climDat_time <- climDat_timeTemp %>% 
  filter(uniqueID %in% goodIDs) %>% 
  dplyr::select(tmin_meanAnnAvg_CLIM:durationFrostFreeDays_meanAnnAvg_3yrAnom, NA_L1CODE, 
         NA_L1NAME, NA_L1KEY, year, newRegion, x, y, soilDepth:totalAvailableWaterHoldingCapacity, uniqueID, spatialID) %>% 
rename("tmin" = tmin_meanAnnAvg_CLIM, 
     "tmax" = tmax_meanAnnAvg_CLIM, #1 
     "tmean" = tmean_meanAnnAvg_CLIM, 
     "prcp" = prcp_meanAnnTotal_CLIM, 
     "t_warm" = T_warmestMonth_meanAnnAvg_CLIM,
     "t_cold" = T_coldestMonth_meanAnnAvg_CLIM, 
     "prcp_wet" = precip_wettestMonth_meanAnnAvg_CLIM,
     "prcp_dry" = precip_driestMonth_meanAnnAvg_CLIM, 
     "prcp_seasonality" = precip_Seasonality_meanAnnAvg_CLIM, #2
     "prcpTempCorr" = PrecipTempCorr_meanAnnAvg_CLIM,  #3
     "abvFreezingMonth" = aboveFreezing_month_meanAnnAvg_CLIM, 
     "isothermality" = isothermality_meanAnnAvg_CLIM, #4
     "annWatDef" = annWaterDeficit_meanAnnAvg_CLIM, 
     "annWetDegDays" = annWetDegDays_meanAnnAvg_CLIM,
     "VPD_mean" = annVPD_mean_meanAnnAvg_CLIM, 
     "VPD_max" = annVPD_max_meanAnnAvg_CLIM, #5
     "VPD_min" = annVPD_min_meanAnnAvg_CLIM, #6
     "VPD_max_95" = annVPD_max_95percentile_CLIM, 
     "annWatDef_95" = annWaterDeficit_95percentile_CLIM, 
     "annWetDegDays_5" = annWetDegDays_5percentile_CLIM, 
     "frostFreeDays_5" = durationFrostFreeDays_5percentile_CLIM, 
     "frostFreeDays" = durationFrostFreeDays_meanAnnAvg_CLIM, 
     "soilDepth" = soilDepth, #7
     "clay" = surfaceClay_perc, 
     "sand" = avgSandPerc_acrossDepth, #8
     "coarse" = avgCoarsePerc_acrossDepth, #9
     "carbon" = avgOrganicCarbonPerc_0_3cm, #10
     "AWHC" = totalAvailableWaterHoldingCapacity,
     ## anomaly variables
     tmean_anom = tmean_meanAnnAvg_3yrAnom, #15
     tmin_anom = tmin_meanAnnAvg_3yrAnom, #16
     tmax_anom = tmax_meanAnnAvg_3yrAnom, #17
    prcp_anom = prcp_meanAnnTotal_3yrAnom, #18
      t_warm_anom = T_warmestMonth_meanAnnAvg_3yrAnom,  #19
     t_cold_anom = T_coldestMonth_meanAnnAvg_3yrAnom, #20
      prcp_wet_anom = precip_wettestMonth_meanAnnAvg_3yrAnom, #21
      precp_dry_anom = precip_driestMonth_meanAnnAvg_3yrAnom,  #22
    prcp_seasonality_anom = precip_Seasonality_meanAnnAvg_3yrAnom, #23 
     prcpTempCorr_anom = PrecipTempCorr_meanAnnAvg_3yrAnom, #24
      aboveFreezingMonth_anom = aboveFreezing_month_meanAnnAvg_3yrAnom, #25  
    isothermality_anom = isothermality_meanAnnAvg_3yrAnom, #26
       annWatDef_anom = annWaterDeficit_meanAnnAvg_3yrAnom, #27
     annWetDegDays_anom = annWetDegDays_meanAnnAvg_3yrAnom,  #28
      VPD_mean_anom = annVPD_mean_meanAnnAvg_3yrAnom, #29
      VPD_min_anom = annVPD_min_meanAnnAvg_3yrAnom,  #30
      VPD_max_anom = annVPD_max_meanAnnAvg_3yrAnom,  #31
     VPD_max_95_anom = annVPD_max_95percentile_3yrAnom, #32
      annWatDef_95_anom = annWaterDeficit_95percentile_3yrAnom, #33 
      annWetDegDays_5_anom = annWetDegDays_5percentile_3yrAnom ,  #34
    frostFreeDays_5_anom = durationFrostFreeDays_5percentile_3yrAnom, #35 
      frostFreeDays_anom = durationFrostFreeDays_meanAnnAvg_3yrAnom #36
  ) %>% 
  dplyr::select(-c(tmin_meanAnnAvg_3yr:Start_3yr)) 


# rename
climDat <- climDat_temp %>% 
  filter(year == 2016) %>% 
  dplyr::select(tmin_meanAnnAvg_CLIM:durationFrostFreeDays_meanAnnAvg_3yrAnom, NA_L1CODE, 
         NA_L1NAME, NA_L1KEY, newRegion, x, y, soilDepth:totalAvailableWaterHoldingCapacity) %>% 
rename("tmin" = tmin_meanAnnAvg_CLIM, 
     "tmax" = tmax_meanAnnAvg_CLIM, #1 
     "tmean" = tmean_meanAnnAvg_CLIM, 
     "prcp" = prcp_meanAnnTotal_CLIM, 
     "t_warm" = T_warmestMonth_meanAnnAvg_CLIM,
     "t_cold" = T_coldestMonth_meanAnnAvg_CLIM, 
     "prcp_wet" = precip_wettestMonth_meanAnnAvg_CLIM,
     "prcp_dry" = precip_driestMonth_meanAnnAvg_CLIM, 
     "prcp_seasonality" = precip_Seasonality_meanAnnAvg_CLIM, #2
     "prcpTempCorr" = PrecipTempCorr_meanAnnAvg_CLIM,  #3
     "abvFreezingMonth" = aboveFreezing_month_meanAnnAvg_CLIM, 
     "isothermality" = isothermality_meanAnnAvg_CLIM, #4
     "annWatDef" = annWaterDeficit_meanAnnAvg_CLIM, 
     "annWetDegDays" = annWetDegDays_meanAnnAvg_CLIM,
     "VPD_mean" = annVPD_mean_meanAnnAvg_CLIM, 
     "VPD_max" = annVPD_max_meanAnnAvg_CLIM, #5
     "VPD_min" = annVPD_min_meanAnnAvg_CLIM, #6
     "VPD_max_95" = annVPD_max_95percentile_CLIM, 
     "annWatDef_95" = annWaterDeficit_95percentile_CLIM, 
     "annWetDegDays_5" = annWetDegDays_5percentile_CLIM, 
     "frostFreeDays_5" = durationFrostFreeDays_5percentile_CLIM, 
     "frostFreeDays" = durationFrostFreeDays_meanAnnAvg_CLIM, 
     "soilDepth" = soilDepth, #7
     "clay" = surfaceClay_perc, 
     "sand" = avgSandPerc_acrossDepth, #8
     "coarse" = avgCoarsePerc_acrossDepth, #9
     "carbon" = avgOrganicCarbonPerc_0_3cm, #10
     "AWHC" = totalAvailableWaterHoldingCapacity,
     ## anomaly variables
     tmean_anom = tmean_meanAnnAvg_3yrAnom, #15
     tmin_anom = tmin_meanAnnAvg_3yrAnom, #16
     tmax_anom = tmax_meanAnnAvg_3yrAnom, #17
    prcp_anom = prcp_meanAnnTotal_3yrAnom, #18
      t_warm_anom = T_warmestMonth_meanAnnAvg_3yrAnom,  #19
     t_cold_anom = T_coldestMonth_meanAnnAvg_3yrAnom, #20
      prcp_wet_anom = precip_wettestMonth_meanAnnAvg_3yrAnom, #21
      precp_dry_anom = precip_driestMonth_meanAnnAvg_3yrAnom,  #22
    prcp_seasonality_anom = precip_Seasonality_meanAnnAvg_3yrAnom, #23 
     prcpTempCorr_anom = PrecipTempCorr_meanAnnAvg_3yrAnom, #24
      aboveFreezingMonth_anom = aboveFreezing_month_meanAnnAvg_3yrAnom, #25  
    isothermality_anom = isothermality_meanAnnAvg_3yrAnom, #26
       annWatDef_anom = annWaterDeficit_meanAnnAvg_3yrAnom, #27
     annWetDegDays_anom = annWetDegDays_meanAnnAvg_3yrAnom,  #28
      VPD_mean_anom = annVPD_mean_meanAnnAvg_3yrAnom, #29
      VPD_min_anom = annVPD_min_meanAnnAvg_3yrAnom,  #30
      VPD_max_anom = annVPD_max_meanAnnAvg_3yrAnom,  #31
     VPD_max_95_anom = annVPD_max_95percentile_3yrAnom, #32
      annWatDef_95_anom = annWaterDeficit_95percentile_3yrAnom, #33 
      annWetDegDays_5_anom = annWetDegDays_5percentile_3yrAnom ,  #34
    frostFreeDays_5_anom = durationFrostFreeDays_5percentile_3yrAnom, #35 
      frostFreeDays_anom = durationFrostFreeDays_meanAnnAvg_3yrAnom #36
  ) %>% 
  dplyr::select(-c(tmin_meanAnnAvg_3yr:Start_3yr))

rm(climDat_temp) 
gc()
##             used   (Mb) gc trigger    (Mb) limit (Mb)   max used    (Mb)
## Ncells   3027807  161.8    4932448   263.5         NA    3877882   207.2
## Vcells 950769060 7253.8 2664304804 20327.1      65536 2664034971 20325.0
## now, scale the contemporary climate data for use in models 
# get the scaling factors 
scaleParams <- modDat_1_s %>% 
  #filter(Year == 2016) %>% 
  dplyr::select(tmin_s:AWHC_s) %>% 
  reframe(across(all_of(names(.)), attributes)) 

# apply the scaling factors to the contemporary climate data 
namesToScale <- climDat %>% 
  dplyr::select(tmin:frostFreeDays, tmean_anom:frostFreeDays_anom, soilDepth:AWHC) %>% 
  names()

climDat_scaled <- map(namesToScale, .f = function(x) {
  x_new <- (climDat[,x] - scaleParams[,paste0(x, "_s")]$`scaled:center`)/scaleParams[,paste0(x, "_s")]$`scaled:scale`
  return(data.frame(x_new))
}) %>% 
  purrr::list_cbind()
names(climDat_scaled) <- paste0(namesToScale, "_s")

climDatPred <- climDat %>% 
  dplyr::select(NA_L1CODE:y) %>% 
  cbind(climDat_scaled)
names(climDatPred)[7:56] <- str_remove(names(climDatPred)[7:56], pattern = "_s$")

# apply the scaling factors to the contemporary climate data used for predictions over time 
climDat_scaled <- map(namesToScale, .f = function(x) {
  x_new <- (climDat_time[,x] - scaleParams[,paste0(x, "_s")]$`scaled:center`)/scaleParams[,paste0(x, "_s")]$`scaled:scale`
  return(data.frame(x_new))
}) %>% 
  purrr::list_cbind()
names(climDat_scaled) <- paste0(namesToScale, "_s")

climDat_time_Pred <- climDat_time %>% 
  dplyr::select(NA_L1CODE:y,spatialID) %>% 
  cbind(climDat_scaled)
names(climDat_time_Pred)[9:58] <- str_remove(names(climDat_time_Pred)[9:58], pattern = "_s$")

rm(climDat_scaled) 
gc()
##             used   (Mb) gc trigger    (Mb) limit (Mb)   max used    (Mb)
## Ncells   3035195  162.1    4932448   263.5         NA    4932448   263.5
## Vcells 976334070 7448.9 2664304804 20327.1      65536 2664034971 20325.0
prednames_s <-  modDat_1_s %>%
  dplyr::select(tmin_s:AWHC_s) %>%
  names()
prednames <- str_replace(prednames_s, pattern = "_s$", replacement = "")

## read in scaled data for future climate 
 forecastClimSoilsDatPred_1 <- readRDS(file = "../../../Data_processed/CoverData/IntermediateAnalysisFiles/Final_ForecastedClimateDataAndSoilsDataForPredictions_BNU-ESM_rcp8_5.rds")
 forecastClimSoilsDatPred_2 <- readRDS(file = "../../../Data_processed/CoverData/IntermediateAnalysisFiles/Final_ForecastedClimateDataAndSoilsDataForPredictions_IPSL-CM5A-MR_rcp8_5.rds")
  # read in unscaled data for future climate
 forecastClimSoilsDat_1 <- readRDS("../../../Data_processed/CoverData/IntermediateAnalysisFiles/Final_ForecastedClimateDataAndSoilsDataForPredictions_BNU-ESM_rcp8_5_UNSCALED.rds")
forecastClimSoilsDat_2 <- readRDS("../../../Data_processed/CoverData/IntermediateAnalysisFiles/Final_ForecastedClimateDataAndSoilsDataForPredictions_IPSL-CM5A-MR_rcp8_5_UNSCALED.rds")

Now, predict using the entire hierarchy of models, starting with the ecoregion classification model

These are predictions of ecoregion classification using contemporary and modeled future climate data

ecoMod <- readRDS("../../ecoregionClassification/ModelFitting/ModelObjects/EcoregionClassificationModel.rds")

## print out model statement
coefficientTable <- data.frame(coefficients(ecoMod))
responseVar <- "P(forest)"
coefficientTable$coefficientName <- rownames(coefficientTable)
coefficientTable$coefficients.ecoMod. <- round(coefficientTable$coefficients.ecoMod., 4)
  
  # print out coefficients w/ coefficient names
tempNames <- paste0(
  apply(coefficientTable, MARGIN = 1, FUN = function(x) {
    if (x["coefficientName"] == "(Intercept)") {
      paste0(x["coefficients.ecoMod."])
    } else {  
      paste0(x["coefficients.ecoMod."], "*", x["coefficientName"])
    }
    }
  ),
  collapse = " + "
)

# print the unscaled model statement
  unscaledModelName <- paste0(responseVar, "~ 1/ (1 + exp(-(", tempNames,")))")
 print(unscaledModelName)
## [1] "P(forest)~ 1/ (1 + exp(-(10.1438 + -0.3121*T_warmestMonth_meanAnnAvg_CLIM +  0.2527*T_coldestMonth_meanAnnAvg_CLIM +  0.0106*precip_wettestMonth_meanAnnAvg_CLIM + -0.0585*annWaterDeficit_meanAnnAvg_CLIM + -2.5967*PrecipTempCorr_meanAnnAvg_CLIM +  0.0582*isothermality_meanAnnAvg_CLIM + -0.0081*soilDepth +  0.0328*avgSandPerc_acrossDepth +  0.0304*avgCoarsePerc_acrossDepth +  0.2734*avgOrganicCarbonPerc_0_3cm)))"

Here are maps of the ecoregion model predictions with both contemporary and modeled future climate data

## get model data used for ecoregion fitting
modDat_ecoregionFit <- readRDS("../../../Data_processed/EcoRegion_climSoilData.rds")

modDat_ecoregionFit <- modDat_ecoregionFit %>% 
  rename("Long" = x, "Lat" = y) %>% 
  dplyr::select(c(newRegion, #swe_meanAnnAvg_CLIM:
                  tmin_meanAnnAvg_CLIM:durationFrostFreeDays_meanAnnAvg_CLIM,
                                         soilDepth                  , surfaceClay_perc          ,                
                                         avgSandPerc_acrossDepth    ,  avgCoarsePerc_acrossDepth,                 
                                         avgOrganicCarbonPerc_0_3cm ,  totalAvailableWaterHoldingCapacity,
                                         Long, Lat)) %>% 
  mutate(newRegion = as.factor(newRegion)) %>% 
  drop_na()

# get climate data from dayMet (d.f. is "climDatPred_unscaled")
climDatPred_unscaled <- climDat
names(climDatPred_unscaled)[c(5, 6, 7, 13, 10, 12, 98, 100, 101, 102)] <- c("T_warmestMonth_meanAnnAvg_CLIM", "T_coldestMonth_meanAnnAvg_CLIM", "precip_wettestMonth_meanAnnAvg_CLIM", "annWaterDeficit_meanAnnAvg_CLIM", "PrecipTempCorr_meanAnnAvg_CLIM", "isothermality_meanAnnAvg_CLIM", "soilDepth", "avgSandPerc_acrossDepth", "avgCoarsePerc_acrossDepth", "avgOrganicCarbonPerc_0_3cm")

#rm(climDat) 
#gc()

# predict with contemporary climate data 
preds_byHand <- climDatPred_unscaled %>% 
  mutate(pred = 1/(1 + exp(-( 9.8726 + -0.2999*T_warmestMonth_meanAnnAvg_CLIM +  0.2456*T_coldestMonth_meanAnnAvg_CLIM +  0.0106*precip_wettestMonth_meanAnnAvg_CLIM + -0.0621*annWaterDeficit_meanAnnAvg_CLIM + -2.7863*PrecipTempCorr_meanAnnAvg_CLIM +  0.0540*isothermality_meanAnnAvg_CLIM + -0.0076*soilDepth +  0.0335*avgSandPerc_acrossDepth +  0.0310*avgCoarsePerc_acrossDepth +  0.2726*avgOrganicCarbonPerc_0_3cm))))
#predictionsModel <- predict(object = ecoMod, type = "response", newdata = climDatPred_unscaled)

# predict with modeled future climate data v1
names(forecastClimSoilsDat_1)[c(8, 9, 10, 16, 13, 15, 48, 50, 51, 52)] <- c("T_warmestMonth_meanAnnAvg_CLIM", "T_coldestMonth_meanAnnAvg_CLIM", "precip_wettestMonth_meanAnnAvg_CLIM", "annWaterDeficit_meanAnnAvg_CLIM", "PrecipTempCorr_meanAnnAvg_CLIM", "isothermality_meanAnnAvg_CLIM", "soilDepth", "avgSandPerc_acrossDepth", "avgCoarsePerc_acrossDepth", "avgOrganicCarbonPerc_0_3cm")
preds_future1_byHand <- forecastClimSoilsDat_1 %>% 
  mutate(pred = 1/(1 + exp(-( 9.8726 + -0.2999*T_warmestMonth_meanAnnAvg_CLIM +  0.2456*T_coldestMonth_meanAnnAvg_CLIM +  0.0106*precip_wettestMonth_meanAnnAvg_CLIM + -0.0621*annWaterDeficit_meanAnnAvg_CLIM + -2.7863*PrecipTempCorr_meanAnnAvg_CLIM +  0.0540*isothermality_meanAnnAvg_CLIM + -0.0076*soilDepth +  0.0335*avgSandPerc_acrossDepth +  0.0310*avgCoarsePerc_acrossDepth +  0.2726*avgOrganicCarbonPerc_0_3cm))))

# predict with modeled future climate data v2
names(forecastClimSoilsDat_2)[c(8, 9, 10, 16, 13, 15, 48, 50, 51, 52)] <- c("T_warmestMonth_meanAnnAvg_CLIM", "T_coldestMonth_meanAnnAvg_CLIM", "precip_wettestMonth_meanAnnAvg_CLIM", "annWaterDeficit_meanAnnAvg_CLIM", "PrecipTempCorr_meanAnnAvg_CLIM", "isothermality_meanAnnAvg_CLIM", "soilDepth", "avgSandPerc_acrossDepth", "avgCoarsePerc_acrossDepth", "avgOrganicCarbonPerc_0_3cm")
preds_future2_byHand <- forecastClimSoilsDat_2 %>% 
  mutate(pred = 1/(1 + exp(-( 9.8726 + -0.2999*T_warmestMonth_meanAnnAvg_CLIM +  0.2456*T_coldestMonth_meanAnnAvg_CLIM +  0.0106*precip_wettestMonth_meanAnnAvg_CLIM + -0.0621*annWaterDeficit_meanAnnAvg_CLIM + -2.7863*PrecipTempCorr_meanAnnAvg_CLIM +  0.0540*isothermality_meanAnnAvg_CLIM + -0.0076*soilDepth +  0.0335*avgSandPerc_acrossDepth +  0.0310*avgCoarsePerc_acrossDepth +  0.2726*avgOrganicCarbonPerc_0_3cm))))


## Plot predictions with contemporary data
# make into a raster
plotObs <- preds_byHand %>% 
         #drop_na(paste(response)) %>% 
  #slice_sample(n = 5e4) %>%
  terra::vect(geom = c("x", "y")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = "pred", 
                   fun = mean, na.rm = TRUE) 

# get the extent of this particular raster, and crop it accordingly

# get the extent of this particular raster, and crop it accordingly
tempExt <- crds(plotObs, na.rm = TRUE)

plotObs_2 <- plotObs %>% 
  crop(ext(min(tempExt[,1]), max(tempExt[,1]),
           min(tempExt[,2]), max(tempExt[,2])) 
       )

map_ecoRegion_current <- ggplot() +
  geom_spatraster(data = plotObs_2) + 
  geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Ecoregion classification predictions using 
                    contemporary dayMet data")) +
  scale_fill_gradient2(low = "brown",
                       mid = "wheat" ,
                       high = "darkgreen" , 
                       midpoint = 0, limits = c(0,1),  na.value = "lightgrey",
                       labels = c("grassShrub", "0.25", "0.50", "0.75", "forest")) + 
  xlim(st_bbox(plotObs_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_2)[c(2,4)])

# for future v1
plotObs_future1 <- preds_future1_byHand %>% 
         #drop_na(paste(response)) %>% 
  #slice_sample(n = 5e4) %>%
  terra::vect(geom = c("x", "y")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = "pred", 
                   fun = mean, na.rm = TRUE)  
    
  
# get the extent of this particular raster, and crop it accordingly

plotObs_future1_2 <- plotObs_future1 %>% 
  crop(ext(min(tempExt[,1]), max(tempExt[,1]),
           min(tempExt[,2]), max(tempExt[,2])) 
       )

map_ecoRegion_future1 <- ggplot() +
  geom_spatraster(data = plotObs_future1_2) + 
  geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_future1_2)),fill=NA ) +
labs(title = paste0("Ecoregion classification predictions using future 
climate data from the BNU-ESM model")) +
  scale_fill_gradient2(low = "brown",
                       mid = "wheat" ,
                       high = "darkgreen" , 
                       midpoint = 0, limits = c(0,1),  na.value = "lightgrey",
                       labels = c("grassShrub", "0.25", "0.50", "0.75", "forest")) + 
  xlim(st_bbox(plotObs_future1_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_future1_2)[c(2,4)])

# for future v2
plotObs_future2 <- preds_future2_byHand %>% 
         #drop_na(paste(response)) %>% 
  #slice_sample(n = 5e4) %>%
  terra::vect(geom = c("x", "y")) %>% 
  terra::set.crs(crs(test_rast)) %>% 
  terra::rasterize(y = test_rast, 
                   field = "pred", 
                   fun = mean, na.rm = TRUE)  
    
  

# get the extent of this particular raster, and crop it accordingly

plotObs_future2_2 <- plotObs_future2 %>% 
  crop(ext(min(tempExt[,1]), max(tempExt[,1]),
           min(tempExt[,2]), max(tempExt[,2])) 
       )

map_ecoRegion_future2 <- ggplot() +
  geom_spatraster(data = plotObs_future2_2) + 
  geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
  geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_future2_2)),fill=NA ) +
labs(title = paste0("Ecoregion classification predictions using future 
climate data from the IPSL-CM5A-MR model")) +
  scale_fill_gradient2(low = "brown",
                       mid = "wheat" ,
                       high = "darkgreen" , 
                       midpoint = 0, limits = c(0,1),  na.value = "lightgrey",
                       labels = c("grassShrub", "0.25", "0.50", "0.75", "forest")) + 
  xlim(st_bbox(plotObs_future2_2)[c(1,3)]) + 
  ylim(st_bbox(plotObs_future2_2)[c(2,4)])

## Plot the maps together
  ggarrange(map_ecoRegion_current, map_ecoRegion_future1, map_ecoRegion_future2, nrow = 1) %>% 
  annotate_figure(fig.lab = "Model Predictions of Ecoregion Classification with Contemporary and Forecasted Climate Data", fig.lab.size = 20)

get model objects for cover models

# Total herbaceous cover: 
# Grass/shrub: 1SE lambda model
totalHerb_GS <-  readRDS("./models/betaLASSO/TotalHerbaceousCover_shrubGrass_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_oneSELambdaGLM.rds")
# Forest: best lambda model
totalHerb_F <-  readRDS("./models/betaLASSO/TotalHerbaceousCover_forest_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")


# Total tree cover: 
# Grass/shrub: best lambda model 
totalTree_GS <- readRDS("./models/betaLASSO/TotalTreeCover_shrubGrass_noTLP_FALSE_removeAnomaliesTRUE_bestLambdaGLM.rds")
# Forest: best lambda model 
totalTree_F <- readRDS("./models/betaLASSO/TotalTreeCover_forest_noTLP_FALSE_removeAnomaliesTRUE_halfSELambdaGLM.rds")
  
# Bare ground: 
# CONUS - 1SE lambda model
bareGround_CONUS <-  readRDS("./models/betaLASSO/BareGroundCover_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_oneSELambdaGLM.rds")

# Shrub: 
# CONUS - best lambda model
shrub_CONUS <- readRDS("./models/betaLASSO/ShrubCover_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")

# Broad-leaved tree: 
# Grass/shrub: best lambda model
BLtree_GS <- readRDS("./models/betaLASSO/AngioTreeCover_prop_shrubGrass_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# Forest:  best lambda model
BLtree_F <- readRDS("./models/betaLASSO/AngioTreeCover_prop_forest_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
  
# Needle-leaved tree: 
# Grass/shrub: best lambda model
NLtree_GS <-  readRDS("./models/betaLASSO/ConifTreeCover_prop_shrubGrass_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# Forest:  best lambda model
NLtree_F <- readRDS("./models/betaLASSO/ConifTreeCover_prop_forest_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")

# Forb:
# CONUS: best lambda model
forb_CONUS <-  readRDS("./models/betaLASSO/ForbCover_prop_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")

# C3 grass:
# CONUS: 1SE lambda model
C3grass_CONUS <- readRDS("./models/betaLASSO/C3GramCover_prop_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_oneSELambdaGLM.rds")

# C4 grass: 
# CONUS: best lambda model
C4grass_CONUS <- readRDS("./models/betaLASSO/C4GramCover_prop_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")

Now, predict absolute cover with scaled climate/weather/soils data

## with contemporary climate data
totalHerb_F_predsContemp <- makePredictions(predictionDF = climDatPred,
                                            modelObject = totalHerb_F) %>% 
  rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsContemp <- makePredictions(predictionDF = climDatPred,
                                            modelObject = totalHerb_GS) %>% 
  rename(absTotalHerb_GS = "modelPreds")

totalTree_F_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = totalTree_F) %>% 
  rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = totalTree_GS) %>% 
  rename(absTotalTree_GS = "modelPreds")

bareGround_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = bareGround_CONUS) %>% 
  rename(absBareGround_CONUS = "modelPreds")

shrub_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = shrub_CONUS) %>% 
  rename(absShrub_CONUS = "modelPreds")

BLtree_F_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = BLtree_F) %>% 
  rename(percBLTree_F = "modelPreds")
BLtree_GS_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = BLtree_GS) %>% 
  rename(percBLTree_GS = "modelPreds")

NLtree_F_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = NLtree_F) %>% 
  rename(percNLTree_F = "modelPreds")
NLtree_GS_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = NLtree_GS) %>% 
  rename(percNLTree_GS = "modelPreds")

C3grass_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = C3grass_CONUS) %>% 
  rename(percC3grass_CONUS = "modelPreds")

C4grass_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = C4grass_CONUS) %>% 
  rename(percC4grass_CONUS = "modelPreds")

forb_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
                                            modelObject = forb_CONUS) %>% 
  rename(percForb_CONUS = "modelPreds")

contempPreds <- totalHerb_F_predsContemp  %>% 
  cbind(totalHerb_GS_predsContemp %>% select(absTotalHerb_GS)) %>% 
  cbind(totalTree_F_predsContemp %>% select(absTotalTree_F)) %>% 
  cbind(totalTree_GS_predsContemp %>% select(absTotalTree_GS)) %>% 
  cbind(bareGround_CONUS_predsContemp %>% select(absBareGround_CONUS)) %>% 
  cbind(shrub_CONUS_predsContemp %>% select(absShrub_CONUS)) %>% 
  cbind(BLtree_F_predsContemp %>% select(percBLTree_F)) %>% 
  cbind(BLtree_GS_predsContemp %>% select(percBLTree_GS)) %>% 
  cbind(NLtree_F_predsContemp %>% select(percNLTree_F)) %>% 
  cbind(NLtree_GS_predsContemp %>% select(percNLTree_GS)) %>% 
  cbind(C3grass_CONUS_predsContemp %>% select(percC3grass_CONUS)) %>% 
  cbind(C4grass_CONUS_predsContemp %>% select(percC4grass_CONUS)) %>% 
  cbind(forb_CONUS_predsContemp %>% select(percForb_CONUS)) %>% 
  cbind(preds_byHand %>% select(pred)) %>% 
  rename(prob_Forest = "pred") %>% 
  mutate(prob_grassShrub = 1 - prob_Forest)

# get predictions with first climate model 
## with contemporary climate data
totalHerb_F_predsFuture1 <- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = totalHerb_F) %>% 
  rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsFuture1 <- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = totalHerb_GS) %>% 
  rename(absTotalHerb_GS = "modelPreds")

totalTree_F_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = totalTree_F) %>% 
  rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = totalTree_GS) %>% 
  rename(absTotalTree_GS = "modelPreds")

bareGround_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = bareGround_CONUS) %>% 
  rename(absBareGround_CONUS = "modelPreds")

shrub_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = shrub_CONUS) %>% 
  rename(absShrub_CONUS = "modelPreds")

BLtree_F_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = BLtree_F) %>% 
  rename(percBLTree_F = "modelPreds")
BLtree_GS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = BLtree_GS) %>% 
  rename(percBLTree_GS = "modelPreds")

NLtree_F_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = NLtree_F) %>% 
  rename(percNLTree_F = "modelPreds")
NLtree_GS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = NLtree_GS) %>% 
  rename(percNLTree_GS = "modelPreds")

C3grass_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = C3grass_CONUS) %>% 
  rename(percC3grass_CONUS = "modelPreds")

C4grass_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = C4grass_CONUS) %>% 
  rename(percC4grass_CONUS = "modelPreds")

forb_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
                                            modelObject = forb_CONUS) %>% 
  rename(percForb_CONUS = "modelPreds")

predsFuture1 <- totalHerb_F_predsFuture1  %>% 
  cbind(totalHerb_GS_predsFuture1 %>% select(absTotalHerb_GS)) %>% 
  cbind(totalTree_F_predsFuture1 %>% select(absTotalTree_F)) %>% 
  cbind(totalTree_GS_predsFuture1 %>% select(absTotalTree_GS)) %>% 
  cbind(bareGround_CONUS_predsFuture1 %>% select(absBareGround_CONUS)) %>% 
  cbind(shrub_CONUS_predsFuture1 %>% select(absShrub_CONUS)) %>% 
  cbind(BLtree_F_predsFuture1 %>% select(percBLTree_F)) %>% 
  cbind(BLtree_GS_predsFuture1 %>% select(percBLTree_GS)) %>% 
  cbind(NLtree_F_predsFuture1 %>% select(percNLTree_F)) %>% 
  cbind(NLtree_GS_predsFuture1 %>% select(percNLTree_GS)) %>% 
  cbind(C3grass_CONUS_predsFuture1 %>% select(percC3grass_CONUS)) %>% 
  cbind(C4grass_CONUS_predsFuture1 %>% select(percC4grass_CONUS)) %>% 
  cbind(forb_CONUS_predsFuture1 %>% select(percForb_CONUS)) %>% 
  cbind(preds_future1_byHand %>% select(pred)) %>% 
  rename(prob_Forest = "pred") %>% 
  mutate(prob_grassShrub = 1 - prob_Forest)

# get predictions with second climate model
totalHerb_F_predsFuture2 <- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = totalHerb_F) %>% 
  rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsFuture2 <- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = totalHerb_GS) %>% 
  rename(absTotalHerb_GS = "modelPreds")

totalTree_F_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = totalTree_F) %>% 
  rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = totalTree_GS) %>% 
  rename(absTotalTree_GS = "modelPreds")

bareGround_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = bareGround_CONUS) %>% 
  rename(absBareGround_CONUS = "modelPreds")

shrub_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = shrub_CONUS) %>% 
  rename(absShrub_CONUS = "modelPreds")

BLtree_F_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = BLtree_F) %>% 
  rename(percBLTree_F = "modelPreds")
BLtree_GS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = BLtree_GS) %>% 
  rename(percBLTree_GS = "modelPreds")

NLtree_F_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = NLtree_F) %>% 
  rename(percNLTree_F = "modelPreds")
NLtree_GS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = NLtree_GS) %>% 
  rename(percNLTree_GS = "modelPreds")

C3grass_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = C3grass_CONUS) %>% 
  rename(percC3grass_CONUS = "modelPreds")

C4grass_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = C4grass_CONUS) %>% 
  rename(percC4grass_CONUS = "modelPreds")

forb_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
                                            modelObject = forb_CONUS) %>% 
  rename(percForb_CONUS = "modelPreds")

predsFuture2 <- totalHerb_F_predsFuture2  %>% 
  cbind(totalHerb_GS_predsFuture2 %>% select(absTotalHerb_GS)) %>% 
  cbind(totalTree_F_predsFuture2 %>% select(absTotalTree_F)) %>% 
  cbind(totalTree_GS_predsFuture2 %>% select(absTotalTree_GS)) %>% 
  cbind(bareGround_CONUS_predsFuture2 %>% select(absBareGround_CONUS)) %>% 
  cbind(shrub_CONUS_predsFuture2 %>% select(absShrub_CONUS)) %>% 
  cbind(BLtree_F_predsFuture2 %>% select(percBLTree_F)) %>% 
  cbind(BLtree_GS_predsFuture2 %>% select(percBLTree_GS)) %>% 
  cbind(NLtree_F_predsFuture2 %>% select(percNLTree_F)) %>% 
  cbind(NLtree_GS_predsFuture2 %>% select(percNLTree_GS)) %>% 
  cbind(C3grass_CONUS_predsFuture2 %>% select(percC3grass_CONUS)) %>% 
  cbind(C4grass_CONUS_predsFuture2 %>% select(percC4grass_CONUS)) %>% 
  cbind(forb_CONUS_predsFuture2 %>% select(percForb_CONUS)) %>% 
  cbind(preds_future2_byHand %>% select(pred)) %>% 
  rename(prob_Forest = "pred") %>% 
  mutate(prob_grassShrub = 1 - prob_Forest)

“Scale” absolute cover predictions according to predicted ecoregion probability

## for contemporary data
contempPreds <- contempPreds %>% 
  mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
         absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F),
         percBLTree_CONUS = (prob_grassShrub * percBLTree_GS) + (prob_Forest * percBLTree_F),
         percNLTree_CONUS = (prob_grassShrub * percNLTree_GS) + (prob_Forest * percNLTree_F)
         )

## for first climate model
predsFuture1 <- predsFuture1 %>% 
  mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
         absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F),
         percBLTree_CONUS = (prob_grassShrub * percBLTree_GS) + (prob_Forest * percBLTree_F),
         percNLTree_CONUS = (prob_grassShrub * percNLTree_GS) + (prob_Forest * percNLTree_F)
         )

## for second climate model 
predsFuture2 <- predsFuture2 %>% 
  mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
         absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F),
         percBLTree_CONUS = (prob_grassShrub * percBLTree_GS) + (prob_Forest * percBLTree_F),
         percNLTree_CONUS = (prob_grassShrub * percNLTree_GS) + (prob_Forest * percNLTree_F)
         )

For percentages (NL/BL tree and C4/C3/Forb), scale so that they sum to one

contempPreds <- contempPreds %>% 
  mutate(percBLTree_scaled_CONUS = percBLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
         percNLTree_scaled_CONUS = percNLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
         percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))

predsFuture1 <- predsFuture1  %>% 
  mutate(percBLTree_scaled_CONUS = percBLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
         percNLTree_scaled_CONUS = percNLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
         percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))

predsFuture2 <- predsFuture2  %>% 
  mutate(percBLTree_scaled_CONUS = percBLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
         percNLTree_scaled_CONUS = percNLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
         percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))

Convert percentages for level 2 cover into absolute cover values

## for contemporary data
contempPreds <- contempPreds %>% 
  mutate(absNLTree_CONUS = (percNLTree_scaled_CONUS * absTotalTree_CONUS), 
         absBLTree_CONUS = (percBLTree_scaled_CONUS * absTotalTree_CONUS), 
         absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS),
         absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS),
         absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS)
  )

## for first climate model
predsFuture1 <- predsFuture1 %>% 
  mutate(absNLTree_CONUS = (percNLTree_scaled_CONUS * absTotalTree_CONUS), 
         absBLTree_CONUS = (percBLTree_scaled_CONUS * absTotalTree_CONUS), 
         absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS),
         absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS),
         absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS)
  )

## for second climate model 
predsFuture2 <- predsFuture2 %>% 
  mutate(absNLTree_CONUS = (percNLTree_scaled_CONUS * absTotalTree_CONUS), 
         absBLTree_CONUS = (percBLTree_scaled_CONUS * absTotalTree_CONUS), 
         absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS),
         absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS),
         absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS)
  )

Show predictions of absolute cover, CONUS-wide (after scaling according to ecoregion classification probability)

# prep observed data for use in figures
obsDat <- modDat_1_s %>% 
  select(Long, Lat, Year, ShrubCover, TotalHerbaceousCover, TotalTreeCover, BareGroundCover, C3GramCover_prop, C4GramCover_prop, ForbCover_prop, AngioTreeCover_prop, ConifTreeCover_prop) %>% 
  mutate(absNLTree_CONUS = ConifTreeCover_prop * TotalTreeCover,
         absBLTree_CONUS = AngioTreeCover_prop * TotalTreeCover, 
         absC3grass_CONUS = C3GramCover_prop * TotalTreeCover, 
         absC4grass_CONUS = C4GramCover_prop * TotalTreeCover,
         absForb_CONUS = ForbCover_prop * TotalTreeCover
         ) %>% 
  rename(absShrub_CONUS = ShrubCover,
         absTotalTree_CONUS = TotalTreeCover, 
         absTotalHerb_CONUS = TotalHerbaceousCover, 
         absBareGround_CONUS = BareGroundCover)

# list of absolute cover names
responseNames <- c("absTotalTree_CONUS", "absTotalHerb_CONUS", "absShrub_CONUS", "absBareGround_CONUS", "absNLTree_CONUS", "absBLTree_CONUS", "absC3grass_CONUS", "absC4grass_CONUS", "absForb_CONUS")

absoluteCoverMaps_contemp <- lapply(responseNames, 
                                      FUN = function(x) {
                                       ## contemporary absolute cover
                                        # rasterize response 
                                        resp_rast <- contempPreds %>% 
                                          terra::vect(geom = c("x", "y")) %>% 
                                          terra::set.crs(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE)   
                                        # get the extent of this particular raster, and crop it accordingly
                                        resp_rast_2 <- resp_rast %>% 
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        # make a map
                                        resp_map <- ggplot() +
                                          geom_spatraster(data = resp_rast_2) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Predictions of absolute cover of ",x,", 
                                                    \n after scaling according to ecoregion prediction")) +
                                          scale_fill_gradient2(low = "brown",
                                                               mid = "wheat" ,
                                                               high = "darkgreen" , 
                                                               midpoint = 0, limits = c(0,1),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist <- ggplot() + 
                                          geom_density(aes(contempPreds[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of absolute cover of ",x)) +
                                          theme_classic()
                                        
                                        ## calculate residuals
                                        # make maps of residuals between absolute cover predictions and observations
                                        # rasterize observations 
                                        plotObservations <- obsDat %>% 
                                          terra::vect(geom = c("Long", "Lat")) %>% 
                                          terra::set.crs(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE) 
                                         plotObservations <- plotObservations %>%  
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        # calculate residuals (observed - predicted)
                                        resids_rast <- plotObservations - resp_rast_2
                                        ## aggregate the residuals raster to make it easier to see 
                                        resids_rastCoarse <- resids_rast %>% 
                                          aggregate(fact = 2, fun = "mean")
                                        # map the residuals
                                        
                                        mapResids <- ggplot() +
                                          geom_spatraster(data = resids_rastCoarse) + 
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% 
                                                    st_crop(st_bbox(resids_rast)),fill=NA ) +
                                          labs(title = paste0("Residuals (observations - predicted absolute cover) for ",x)) +
                                          scale_fill_gradient2(low = "red",
                                                               mid = "white" ,
                                                               high = "blue" , 
                                                               midpoint = 0,   na.value = "lightgrey",
                                                               limits = c(-1,1)
                                          )  + 
                                          xlim(st_bbox(resids_rast)[c(1,3)]) + 
                                          ylim(st_bbox(resids_rast)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram of residuals
                                        histResids <- ggplot() + 
                                          geom_density(aes(terra::values(resids_rast$mean, na.rm = TRUE, mat = FALSE)), 
                                                       col = "black", fill = "darkgrey") + 
                                          geom_vline(aes(xintercept = 0), linetype = 2) +
                                          xlab(paste0("Residuals (obs - pred) ",x)) +
                                          theme_classic()
                                        
                                         ## future model 1 absolute cover
                                        # rasterize response 
                                        resp_rast_future1 <- predsFuture1 %>% 
                                          terra::vect(geom = c("x", "y")) %>% 
                                          terra::set.crs(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE)   
                                        # get the extent of this particular raster, and crop it accordingly
                                        resp_rast_future1_2 <- resp_rast_future1 %>% 
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        # calculate differences between contemporary preds and future 1 preds (deltas)
                                        rast_deltas_future1 <- resp_rast_future1_2 - resp_rast_2
                                        # make a map
                                        resp_map_future1 <- ggplot() +
                                          geom_spatraster(data = rast_deltas_future1) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Future predictions of absolute cover of ",x,", \n compared to contemporary predictions (future - contemprary)  \n - climate model: BNU-ESM")) +
                                          scale_fill_gradient2(low = "orange",
                                                               mid = "white" ,
                                                               high = "purple" , 
                                                               midpoint = 0, limits = c(-1,1),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist_future1 <- ggplot() + 
                                          geom_density(aes(predsFuture1[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of absolute cover of ",x, 
                                                      " \n in future scenario (model = BNU-ESM)")) +
                                          theme_classic()
                                        
                                        ## future model 2 absolute cover
                                        # rasterize response 
                                        resp_rast_future2 <- predsFuture2 %>% 
                                          terra::vect(geom = c("x", "y")) %>% 
                                          terra::set.crs(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE)   
                                        # get the extent of this particular raster, and crop it accordingly
                                        resp_rast_future2_2 <- resp_rast_future2 %>% 
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        # calculate differences between contemporary preds and future 1 preds (deltas)
                                        rast_deltas_future2 <- resp_rast_future2_2 - resp_rast_2
                                        # make a map
                                        resp_map_future2 <- ggplot() +
                                          geom_spatraster(data = rast_deltas_future2) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Future predictions of absolute cover of ",x,", \n compared to contemporary predictions (future - contemprary) \n- climate model: IPSL-CM5A-MR (France)")) +
                                          scale_fill_gradient2(low = "orange",
                                                               mid = "white" ,
                                                               high = "purple" , 
                                                               midpoint = 0, limits = c(-1,1),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist_future2 <- ggplot() + 
                                          geom_density(aes(predsFuture2[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of absolute cover of ",x, 
                                                      " \n in future scenario (model = IPSL-CM5A-MR (France))")) +
                                          theme_classic()
                                        
                                        
                                        return(list("rast" = resp_rast_2, 
                                                    "map" = resp_map, 
                                                    "hist" = resp_hist,
                                                    "rast_resids" = resids_rast,
                                                    "map_resids" = mapResids,
                                                    "hist_resids" = histResids,
                                                    "rast_future1" = resp_rast_future1_2,
                                                    "map_future1" = resp_map_future1, 
                                                    "hist_future1" = resp_hist_future1,
                                                    "rast_future2" = resp_rast_future2_2,
                                                    "map_future2" = resp_map_future2, 
                                                    "hist_future2" = resp_hist_future2))
                                      })


names(absoluteCoverMaps_contemp) <- c("absTotalTree_CONUS", "absTotalHerb_CONUS", "absShrub_CONUS", "absBareGround_CONUS", "absNLTree_CONUS", "absBLTree_CONUS", "absC3grass_CONUS", "absC4grass_CONUS", "absForb_CONUS")
ggarrange(
ggarrange(absoluteCoverMaps_contemp$absTotalTree_CONUS$map, 
          absoluteCoverMaps_contemp$absTotalTree_CONUS$map_resids,
          absoluteCoverMaps_contemp$absTotalTree_CONUS$map_future1,
          absoluteCoverMaps_contemp$absTotalTree_CONUS$map_future2,
          absoluteCoverMaps_contemp$absTotalTree_CONUS$hist,
          absoluteCoverMaps_contemp$absTotalTree_CONUS$hist_resids,
           absoluteCoverMaps_contemp$absTotalTree_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absTotalTree_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(absoluteCoverMaps_contemp$absTotalHerb_CONUS$map, 
          absoluteCoverMaps_contemp$absTotalHerb_CONUS$map_resids,
          absoluteCoverMaps_contemp$absTotalHerb_CONUS$map_future1,
          absoluteCoverMaps_contemp$absTotalHerb_CONUS$map_future2,
          absoluteCoverMaps_contemp$absTotalHerb_CONUS$hist,
          absoluteCoverMaps_contemp$absTotalHerb_CONUS$hist_resids, 
          absoluteCoverMaps_contemp$absTotalHerb_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absTotalHerb_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(absoluteCoverMaps_contemp$absShrub_CONUS$map, 
          absoluteCoverMaps_contemp$absShrub_CONUS$map_resids,
          absoluteCoverMaps_contemp$absShrub_CONUS$map_future1,
          absoluteCoverMaps_contemp$absShrub_CONUS$map_future2,
          absoluteCoverMaps_contemp$absShrub_CONUS$hist,
          absoluteCoverMaps_contemp$absShrub_CONUS$hist_resids,
          absoluteCoverMaps_contemp$absShrub_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absShrub_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(absoluteCoverMaps_contemp$absBareGround_CONUS$map, 
          absoluteCoverMaps_contemp$absBareGround_CONUS$map_resids,
          absoluteCoverMaps_contemp$absBareGround_CONUS$map_future1,
          absoluteCoverMaps_contemp$absBareGround_CONUS$map_future2,
          absoluteCoverMaps_contemp$absBareGround_CONUS$hist,
          absoluteCoverMaps_contemp$absBareGround_CONUS$hist_resids, 
          absoluteCoverMaps_contemp$absBareGround_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absBareGround_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(absoluteCoverMaps_contemp$absNLTree_CONUS$map, 
          absoluteCoverMaps_contemp$absNLTree_CONUS$map_resids,
          absoluteCoverMaps_contemp$absNLTree_CONUS$map_future1,
          absoluteCoverMaps_contemp$absNLTree_CONUS$map_future2,
          absoluteCoverMaps_contemp$absNLTree_CONUS$hist,
          absoluteCoverMaps_contemp$absNLTree_CONUS$hist_resids, 
          absoluteCoverMaps_contemp$absNLTree_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absNLTree_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(absoluteCoverMaps_contemp$absBLTree_CONUS$map, 
          absoluteCoverMaps_contemp$absBLTree_CONUS$map_resids,
          absoluteCoverMaps_contemp$absBLTree_CONUS$map_future1,
          absoluteCoverMaps_contemp$absBLTree_CONUS$map_future2,
          absoluteCoverMaps_contemp$absBLTree_CONUS$hist,
          absoluteCoverMaps_contemp$absBLTree_CONUS$hist_resids,
          absoluteCoverMaps_contemp$absBLTree_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absBLTree_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(absoluteCoverMaps_contemp$absC3grass_CONUS$map, 
          absoluteCoverMaps_contemp$absC3grass_CONUS$map_resids,
          absoluteCoverMaps_contemp$absC3grass_CONUS$map_future1,
          absoluteCoverMaps_contemp$absC3grass_CONUS$map_future2,
          absoluteCoverMaps_contemp$absC3grass_CONUS$hist,
          absoluteCoverMaps_contemp$absC3grass_CONUS$hist_resids, 
          absoluteCoverMaps_contemp$absC3grass_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absC3grass_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(absoluteCoverMaps_contemp$absC4grass_CONUS$map, 
          absoluteCoverMaps_contemp$absC4grass_CONUS$map_resids,
          absoluteCoverMaps_contemp$absC4grass_CONUS$map_future1,
          absoluteCoverMaps_contemp$absC4grass_CONUS$map_future2,
          absoluteCoverMaps_contemp$absC4grass_CONUS$hist,
          absoluteCoverMaps_contemp$absC4grass_CONUS$hist_resids, 
          absoluteCoverMaps_contemp$absC4grass_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absC4grass_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(absoluteCoverMaps_contemp$absForb_CONUS$map, 
          absoluteCoverMaps_contemp$absForb_CONUS$map_resids,
          absoluteCoverMaps_contemp$absForb_CONUS$map_future1,
          absoluteCoverMaps_contemp$absForb_CONUS$map_future2,
          absoluteCoverMaps_contemp$absForb_CONUS$hist,
          absoluteCoverMaps_contemp$absForb_CONUS$hist_resids, 
          absoluteCoverMaps_contemp$absForb_CONUS$hist_future1, 
          absoluteCoverMaps_contemp$absForb_CONUS$hist_future2,
          ncol = 4,
          nrow = 2,
          heights = c(1,.35)
          ), 
ncol = 1
)

Get predictions of absolute cover across time using observed climate/weather/soils data

# get climate data from dayMet (d.f. is "climDatPred_unscaled")
climDat_time_unscaled <- climDat_time
names(climDat_time_unscaled)[c(5, 6, 7, 13, 10, 12, 99, 101, 102, 103)] <- c("T_warmestMonth_meanAnnAvg_CLIM", "T_coldestMonth_meanAnnAvg_CLIM", "precip_wettestMonth_meanAnnAvg_CLIM", "annWaterDeficit_meanAnnAvg_CLIM", "PrecipTempCorr_meanAnnAvg_CLIM", "isothermality_meanAnnAvg_CLIM", "soilDepth", "avgSandPerc_acrossDepth", "avgCoarsePerc_acrossDepth", "avgOrganicCarbonPerc_0_3cm")

# predict with contemporary climate data 
preds_byHand_time <- climDat_time_unscaled %>% 
  mutate(pred = 1/(1 + exp(-( 9.8726 + -0.2999*T_warmestMonth_meanAnnAvg_CLIM +  0.2456*T_coldestMonth_meanAnnAvg_CLIM +  0.0106*precip_wettestMonth_meanAnnAvg_CLIM + -0.0621*annWaterDeficit_meanAnnAvg_CLIM + -2.7863*PrecipTempCorr_meanAnnAvg_CLIM +  0.0540*isothermality_meanAnnAvg_CLIM + -0.0076*soilDepth +  0.0335*avgSandPerc_acrossDepth +  0.0310*avgCoarsePerc_acrossDepth +  0.2726*avgOrganicCarbonPerc_0_3cm))))
#predictionsModel <- predict(object = ecoMod, type = "response", newdata = climDatPred_unscaled)

## with contemporary climate data
totalHerb_F_predsOverTime <- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = totalHerb_F) %>% 
  rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsOverTime <- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = totalHerb_GS) %>% 
  rename(absTotalHerb_GS = "modelPreds")

totalTree_F_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = totalTree_F) %>% 
  rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = totalTree_GS) %>% 
  rename(absTotalTree_GS = "modelPreds")

bareGround_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = bareGround_CONUS) %>% 
  rename(absBareGround_CONUS = "modelPreds")

shrub_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = shrub_CONUS) %>% 
  rename(absShrub_CONUS = "modelPreds")

BLtree_F_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = BLtree_F) %>% 
  rename(percBLTree_F = "modelPreds")
BLtree_GS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = BLtree_GS) %>% 
  rename(percBLTree_GS = "modelPreds")

NLtree_F_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = NLtree_F) %>% 
  rename(percNLTree_F = "modelPreds")
NLtree_GS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = NLtree_GS) %>% 
  rename(percNLTree_GS = "modelPreds")

C3grass_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = C3grass_CONUS) %>% 
  rename(percC3grass_CONUS = "modelPreds")

C4grass_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = C4grass_CONUS) %>% 
  rename(percC4grass_CONUS = "modelPreds")

forb_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
                                            modelObject = forb_CONUS) %>% 
  rename(percForb_CONUS = "modelPreds")

predsOverTime <- totalHerb_F_predsOverTime  %>% 
  cbind(totalHerb_GS_predsOverTime %>% select(absTotalHerb_GS)) %>% 
  cbind(totalTree_F_predsOverTime %>% select(absTotalTree_F)) %>% 
  cbind(totalTree_GS_predsOverTime %>% select(absTotalTree_GS)) %>% 
  cbind(bareGround_CONUS_predsOverTime %>% select(absBareGround_CONUS)) %>% 
  cbind(shrub_CONUS_predsOverTime %>% select(absShrub_CONUS)) %>% 
  cbind(BLtree_F_predsOverTime %>% select(percBLTree_F)) %>% 
  cbind(BLtree_GS_predsOverTime %>% select(percBLTree_GS)) %>% 
  cbind(NLtree_F_predsOverTime %>% select(percNLTree_F)) %>% 
  cbind(NLtree_GS_predsOverTime %>% select(percNLTree_GS)) %>% 
  cbind(C3grass_CONUS_predsOverTime %>% select(percC3grass_CONUS)) %>% 
  cbind(C4grass_CONUS_predsOverTime %>% select(percC4grass_CONUS)) %>% 
  cbind(forb_CONUS_predsOverTime %>% select(percForb_CONUS)) %>% 
  cbind(preds_byHand_time %>% select(pred)) %>% 
  rename(prob_Forest = "pred") %>% 
  mutate(prob_grassShrub = 1 - prob_Forest)

## now, scale predictions according to ecoregion percentage
predsOverTime <- predsOverTime %>% 
  mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
         absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F),
         percBLTree_CONUS = (prob_grassShrub * percBLTree_GS) + (prob_Forest * percBLTree_F),
         percNLTree_CONUS = (prob_grassShrub * percNLTree_GS) + (prob_Forest * percNLTree_F)
         )
## then, scale the percentages so they sum to 1
predsOverTime <- predsOverTime %>% 
  mutate(percBLTree_scaled_CONUS = percBLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
         percNLTree_scaled_CONUS = percNLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
         percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
         percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))

## convert percentages into absolute cover for level 2 cover classes
predsOverTime <- predsOverTime %>% 
    mutate(absNLTree_CONUS = (percNLTree_scaled_CONUS * absTotalTree_CONUS), 
         absBLTree_CONUS = (percBLTree_scaled_CONUS * absTotalTree_CONUS), 
         absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS),
         absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS),
         absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS)
  )

Now, visualize predictions of absolute cover over time

# change spatialIDs to discrete numbers, since the long form is causing some sort of floating point issue? 
spatID_lu <- data.frame("spatialID" = unique(predsOverTime$spatialID), 
                        "uniqueID" = 1:length(predsOverTime$spatialID))
predsOverTime <- predsOverTime %>% 
  left_join(spatID_lu)

# make into a long format
predsOverTime_long <- predsOverTime %>% 
  select(year,x, y, spatialID, uniqueID, absShrub_CONUS, absBareGround_CONUS, absTotalHerb_CONUS, absTotalTree_CONUS, absNLTree_CONUS:absForb_CONUS) %>% 
  pivot_longer(cols = c(absShrub_CONUS, absBareGround_CONUS, absTotalHerb_CONUS, absTotalTree_CONUS, absNLTree_CONUS:absForb_CONUS)) %>% 
  mutate(year = as.integer(year))

# calculate the coefficient of variation for each cover type at each location
predsOverTime_cv <- predsOverTime_long %>% 
  group_by(uniqueID, name) %>% 
  summarize(cv = sd(value)/mean(value))
## add x,y 
predsOverTime_cv <- predsOverTime_cv %>% 
  left_join(predsOverTime %>% select(uniqueID, x, y))

# visualize change over time for a selection of points
predsOverTime_long %>% 
  filter(uniqueID %in% 1:50) %>% 
  filter(name %in% c("absShrub_CONUS", "absBareGround_CONUS",  "absNLTree_CONUS", 
                     "absBLTree_CONUS", "absC3grass_CONUS", "absC4grass_CONUS", "absForb_CONUS")) %>% 
  ggplot() + 
  facet_wrap(~uniqueID) + 
  geom_area(aes(x = year, y = value, fill = name)) + 
  geom_line(data = predsOverTime_long %>% 
              filter(uniqueID %in% 1:50) %>% 
              filter(name %in% c("absTotalHerb_CONUS", "absTotalTree_CONUS")),
            aes(x = year, y = value, linetype = name), col = "black") +
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 90)) + 
  labs(title = "Change in modeled absolute cover by functional type from \n 2010 to 2020 for a subset of locations")

predsOverTime_cv %>%  
  ggplot() + 
  geom_boxplot(aes(x = name, y = cv, col = name)) +
  theme_minimal() + 
  labs(title = "Coefficients of variation of model-predicted absolute \n cover at a given location for each functional group from 2010-2020",
        x = "functional group",
       y = "coefficient of variation") +
  theme(axis.text.x = element_text(angle = 90))

# show the locations of each point that is being used for this temporal analysis
#cropped_states <- cropped_states %>% 
 # rename(geometry = Shape)
 
ggplot() + 
  facet_wrap(~name) + 
  geom_sf(data = cropped_states, aes()) + 
  stat_summary2d(data = predsOverTime_cv, aes(x = x, y = y, z = cv), fun = ~ mean(.x, na.rm = TRUE)) +
  geom_hex(data = predsOverTime_cv, aes(x = x, y = y, col = cv), stat = "mean", na.rm = TRUE) + 
  guides(fill = guide_legend(title = "CV")) + 
  theme_minimal() + 
  labs(title = "Locations of points included in the random \n sample for examining temporal variation \n color-coded by the coefficient of variation ")

Now, relativize the predictions across the different functional groups

First option: relativize cover so that cover for all groups sums to 1

## for contemporary data
contempPreds <- contempPreds %>% 
  mutate(totalAbsCover_CONUS = (absTotalHerb_CONUS + absTotalTree_CONUS + absBareGround_CONUS + absShrub_CONUS)) %>% 
  mutate(relCoverA_totalTree = absTotalTree_CONUS/totalAbsCover_CONUS,
         relCoverA_totalHerb = absTotalHerb_CONUS/totalAbsCover_CONUS,
         relCoverA_shrub = absShrub_CONUS/totalAbsCover_CONUS,
         relCoverA_bareGround = absBareGround_CONUS/totalAbsCover_CONUS,
         relCoverA_NLTree = absNLTree_CONUS/totalAbsCover_CONUS,
         relCoverA_BLTree = absBLTree_CONUS/totalAbsCover_CONUS,
         relCoverA_C3Grass = absC3grass_CONUS/totalAbsCover_CONUS,
         relCoverA_C4Grass = absC4grass_CONUS/totalAbsCover_CONUS,
         relCoverA_Forb = absForb_CONUS/totalAbsCover_CONUS
         )

## for first climate model
predsFuture1 <- predsFuture1 %>% 
  mutate(totalAbsCover_CONUS = (absTotalHerb_CONUS + absTotalTree_CONUS + absBareGround_CONUS + absShrub_CONUS)) %>% 
  mutate(relCoverA_totalTree = absTotalTree_CONUS/totalAbsCover_CONUS,
         relCoverA_totalHerb = absTotalHerb_CONUS/totalAbsCover_CONUS,
         relCoverA_shrub = absShrub_CONUS/totalAbsCover_CONUS,
         relCoverA_bareGround = absBareGround_CONUS/totalAbsCover_CONUS,
         relCoverA_NLTree = absNLTree_CONUS/totalAbsCover_CONUS,
         relCoverA_BLTree = absBLTree_CONUS/totalAbsCover_CONUS,
         relCoverA_C3Grass = absC3grass_CONUS/totalAbsCover_CONUS,
         relCoverA_C4Grass = absC4grass_CONUS/totalAbsCover_CONUS,
         relCoverA_Forb = absForb_CONUS/totalAbsCover_CONUS
         )

## for second climate model 
predsFuture2 <- predsFuture2 %>% 
  mutate(totalAbsCover_CONUS = (absTotalHerb_CONUS + absTotalTree_CONUS + absBareGround_CONUS + absShrub_CONUS)) %>% 
  mutate(relCoverA_totalTree = absTotalTree_CONUS/totalAbsCover_CONUS,
         relCoverA_totalHerb = absTotalHerb_CONUS/totalAbsCover_CONUS,
         relCoverA_shrub = absShrub_CONUS/totalAbsCover_CONUS,
         relCoverA_bareGround = absBareGround_CONUS/totalAbsCover_CONUS,
         relCoverA_NLTree = absNLTree_CONUS/totalAbsCover_CONUS,
         relCoverA_BLTree = absBLTree_CONUS/totalAbsCover_CONUS,
         relCoverA_C3Grass = absC3grass_CONUS/totalAbsCover_CONUS,
         relCoverA_C4Grass = absC4grass_CONUS/totalAbsCover_CONUS,
         relCoverA_Forb = absForb_CONUS/totalAbsCover_CONUS
         )

Now, show the results in map form

# list of absolute cover names
responseNames <- c("relCoverA_totalTree", "relCoverA_totalHerb", "relCoverA_shrub", "relCoverA_bareGround",
                   "relCoverA_NLTree", "relCoverA_BLTree", "relCoverA_C3Grass", "relCoverA_C4Grass", "relCoverA_Forb",
                   "totalAbsCover_CONUS")

relCover_A_Maps_contemp <- lapply(responseNames, 
                                      FUN = function(x) {
                                       ## contemporary absolute cover
                                        # rasterize response 
                                        resp_rast <- contempPreds %>% 
                                          terra::vect(geom = c("x", "y")) %>% 
                                          terra::set.crs(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE)   
                                        # get the extent of this particular raster, and crop it accordingly
                                        resp_rast_2 <- resp_rast %>% 
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        
                                        
                                         ## future model 1 absolute cover
                                        # rasterize response 
                                        resp_rast_future1 <- predsFuture1 %>% 
                                          terra::vect(geom = c("x", "y")) %>% 
                                          terra::set.crs(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE)   
                                        # get the extent of this particular raster, and crop it accordingly
                                        resp_rast_future1_2 <- resp_rast_future1 %>% 
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        # calculate differences between contemporary preds and future 1 preds (deltas)
                                        rast_deltas_future1 <- resp_rast_future1_2 - resp_rast_2
                                        
                                        ## future model 2 absolute cover
                                        # rasterize response 
                                        resp_rast_future2 <- predsFuture2 %>% 
                                          terra::vect(geom = c("x", "y")) %>% 
                                          terra::set.crs(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE)   
                                        # get the extent of this particular raster, and crop it accordingly
                                        resp_rast_future2_2 <- resp_rast_future2 %>% 
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        # calculate differences between contemporary preds and future 1 preds (deltas)
                                        rast_deltas_future2 <- resp_rast_future2_2 - resp_rast_2
                                        
                                        
                                        if (x == "totalAbsCover_CONUS") {
                                          # make a map of predictions w/ current data 
                                        resp_map <- ggplot() +
                                          geom_spatraster(data = resp_rast_2) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Predictions of Total Absolute Cover with contemporary data")) +
                                          scale_fill_gradient2(low = "lightblue",
                                                               #mid = "wheat" ,
                                                               high = "darkblue" , 
                                                               #midpoint = 0, 
                                                               limits = c(0,max(terra::values(resp_rast_2$mean), na.rm= TRUE)),  
                                                              na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist <- ggplot() + 
                                          geom_density(aes(contempPreds[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of Total Absolute Cover with contemporary data")) +
                                          theme_classic()
                                        
                                         # make a map of future preds w/ model 1
                                         
                                        resp_map_future1 <- ggplot() +
                                          geom_spatraster(data = rast_deltas_future1) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Future predictions of Total Absolute Cover \n compared to contemporary predictions (future - contemprary)  \n - climate model: BNU-ESM")) +
                                          scale_fill_gradient2(low = "orange",
                                                               mid = "white" ,
                                                               high = "purple" , 
                                                               midpoint = 0, limits = c(-2,2),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist_future1 <- ggplot() + 
                                          geom_density(aes(predsFuture1[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of relative cover of ",x, 
                                                      " \n in future scenario (model = BNU-ESM)")) +
                                          theme_classic()
                                        
                                        # make a map of future preds w/ model 2
                                        resp_map_future2 <- ggplot() +
                                          geom_spatraster(data = rast_deltas_future2) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Future predictions of Total Absolute Cover, \n compared to contemporary predictions (future - contemprary) \n- climate model: IPSL-CM5A-MR (France)")) +
                                          scale_fill_gradient2(low = "orange",
                                                               mid = "white" ,
                                                               high = "purple" , 
                                                               midpoint = 0, limits = c(-2,2),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist_future2 <- ggplot() + 
                                          geom_density(aes(predsFuture2[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of relative cover of ",x, 
                                                      " \n in future scenario (model = IPSL-CM5A-MR (France))")) +
                                          theme_classic()
                                        
                                        } else {
                                          # make a map of predictions w/ current data 
                                        resp_map <- ggplot() +
                                          geom_spatraster(data = resp_rast_2) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Predictions of relative cover of ",x,", 
                                                    \n when all functional groups sum to 1")) +
                                          scale_fill_gradient2(low = "brown",
                                                               mid = "wheat" ,
                                                               high = "darkgreen" , 
                                                               midpoint = 0, limits = c(0,1),  
                                                              na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist <- ggplot() + 
                                          geom_density(aes(contempPreds[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of relative cover of ",x)) +
                                          theme_classic()
                                        
                                         # make a map of future preds w/ model 1
                                         
                                        resp_map_future1 <- ggplot() +
                                          geom_spatraster(data = rast_deltas_future1) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Future predictions of relative cover of ",x,", \n compared to contemporary predictions (future - contemprary)  \n - climate model: BNU-ESM")) +
                                          scale_fill_gradient2(low = "orange",
                                                               mid = "white" ,
                                                               high = "purple" , 
                                                               midpoint = 0, limits = c(-1,1),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist_future1 <- ggplot() + 
                                          geom_density(aes(predsFuture1[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of relative cover of ",x, 
                                                      " \n in future scenario (model = BNU-ESM)")) +
                                          theme_classic()
                                        
                                        # make a map of future preds w/ model 2
                                        resp_map_future2 <- ggplot() +
                                          geom_spatraster(data = rast_deltas_future2) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Future predictions of relative cover of ",
                                                                        x,
                                                                        ", \n compared to contemporary predictions (future - contemprary) \n- climate model: IPSL-CM5A-MR (France)")) +
                                          scale_fill_gradient2(low = "orange",
                                                               mid = "white" ,
                                                               high = "purple" , 
                                                               midpoint = 0, limits = c(-1,1),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist_future2 <- ggplot() + 
                                          geom_density(aes(predsFuture2[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of relative cover of ",x, 
                                                      " \n in future scenario (model = IPSL-CM5A-MR (France))")) +
                                          theme_classic()
                                        }
                                        
                                        
                                        return(list("rast" = resp_rast_2, 
                                                    "map" = resp_map, 
                                                    "hist" = resp_hist,
                                                    "rast_future1" = resp_rast_future1_2,
                                                    "map_future1" = resp_map_future1, 
                                                    "hist_future1" = resp_hist_future1,
                                                    "rast_future2" = resp_rast_future2_2,
                                                    "map_future2" = resp_map_future2, 
                                                    "hist_future2" = resp_hist_future2))
                                      })


names(relCover_A_Maps_contemp) <- c("relCoverA_totalTree", "relCoverA_totalHerb", "relCoverA_shrub", "relCoverA_bareGround",
                   "relCoverA_NLTree", "relCoverA_BLTree", "relCoverA_C3Grass", "relCoverA_C4Grass", "relCoverA_Forb",
                   "totalAbsCover_CONUS")
ggarrange(
ggarrange(relCover_A_Maps_contemp$relCoverA_totalTree$map, 
          relCover_A_Maps_contemp$relCoverA_totalTree$map_future1,
          relCover_A_Maps_contemp$relCoverA_totalTree$map_future2,
          relCover_A_Maps_contemp$relCoverA_totalTree$hist,
           relCover_A_Maps_contemp$relCoverA_totalTree$hist_future1, relCover_A_Maps_contemp$relCoverA_totalTree$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_A_Maps_contemp$relCoverA_totalHerb$map, 
          relCover_A_Maps_contemp$relCoverA_totalHerb$map_future1,
          relCover_A_Maps_contemp$relCoverA_totalHerb$map_future2,
          relCover_A_Maps_contemp$relCoverA_totalHerb$hist,
          relCover_A_Maps_contemp$relCoverA_totalHerb$hist_future1, 
          relCover_A_Maps_contemp$relCoverA_totalHerb$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_A_Maps_contemp$relCoverA_shrub$map, 
          relCover_A_Maps_contemp$relCoverA_shrub$map_future1,
          relCover_A_Maps_contemp$relCoverA_shrub$map_future2,
          relCover_A_Maps_contemp$relCoverA_shrub$hist,
          relCover_A_Maps_contemp$relCoverA_shrub$hist_future1, 
          relCover_A_Maps_contemp$relCoverA_shrub$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_A_Maps_contemp$relCoverA_bareGround$map, 
          relCover_A_Maps_contemp$relCoverA_bareGround$map_future1,
          relCover_A_Maps_contemp$relCoverA_bareGround$map_future2,
          relCover_A_Maps_contemp$relCoverA_bareGround$hist,
          relCover_A_Maps_contemp$relCoverA_bareGround$hist_future1, 
          relCover_A_Maps_contemp$relCoverA_bareGround$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_A_Maps_contemp$relCoverA_NLTree$map, 
          relCover_A_Maps_contemp$relCoverA_NLTree$map_future1,
          relCover_A_Maps_contemp$relCoverA_NLTree$map_future2,
          relCover_A_Maps_contemp$relCoverA_NLTree$hist,
          relCover_A_Maps_contemp$relCoverA_NLTree$hist_future1, 
          relCover_A_Maps_contemp$relCoverA_NLTree$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_A_Maps_contemp$relCoverA_BLTree$map, 
          relCover_A_Maps_contemp$relCoverA_BLTree$map_future1,
          relCover_A_Maps_contemp$relCoverA_BLTree$map_future2,
          relCover_A_Maps_contemp$relCoverA_BLTree$hist,
          relCover_A_Maps_contemp$relCoverA_BLTree$hist_future1, 
          relCover_A_Maps_contemp$relCoverA_BLTree$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_A_Maps_contemp$relCoverA_C3Grass$map, 
          relCover_A_Maps_contemp$relCoverA_C3Grass$map_future1,
          relCover_A_Maps_contemp$relCoverA_C3Grass$map_future2,
          relCover_A_Maps_contemp$relCoverA_C3Grass$hist,
          relCover_A_Maps_contemp$relCoverA_C3Grass$hist_future1, 
          relCover_A_Maps_contemp$relCoverA_C3Grass$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_A_Maps_contemp$relCoverA_C4Grass$map, 
          relCover_A_Maps_contemp$relCoverA_C4Grass$map_future1,
          relCover_A_Maps_contemp$relCoverA_C4Grass$map_future2,
          relCover_A_Maps_contemp$relCoverA_C4Grass$hist,
          relCover_A_Maps_contemp$relCoverA_C4Grass$hist_future1, 
          relCover_A_Maps_contemp$relCoverA_C4Grass$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_A_Maps_contemp$relCoverA_Forb$map, 
          relCover_A_Maps_contemp$relCoverA_Forb$map_future1,
          relCover_A_Maps_contemp$relCoverA_Forb$map_future2,
          relCover_A_Maps_contemp$relCoverA_Forb$hist,
          relCover_A_Maps_contemp$relCoverA_Forb$hist_future1, 
          relCover_A_Maps_contemp$relCoverA_Forb$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ), 
ggarrange(relCover_A_Maps_contemp$totalAbsCover_CONUS$map, 
          relCover_A_Maps_contemp$totalAbsCover_CONUS$map_future1,
          relCover_A_Maps_contemp$totalAbsCover_CONUS$map_future2,
          relCover_A_Maps_contemp$totalAbsCover_CONUS$hist,
          relCover_A_Maps_contemp$totalAbsCover_CONUS$hist_future1, 
          relCover_A_Maps_contemp$totalAbsCover_CONUS$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ncol = 1
)

Now, show how the relativized cover values change over time

# calculate relativized cover (version w/ all functional groups summing to 1)
predsOverTime <- predsOverTime %>% 
mutate(totalAbsCover_CONUS = (absTotalHerb_CONUS + absTotalTree_CONUS + absBareGround_CONUS + absShrub_CONUS)) %>% 
  mutate(relCoverA_totalTree = absTotalTree_CONUS/totalAbsCover_CONUS,
         relCoverA_totalHerb = absTotalHerb_CONUS/totalAbsCover_CONUS,
         relCoverA_shrub = absShrub_CONUS/totalAbsCover_CONUS,
         relCoverA_bareGround = absBareGround_CONUS/totalAbsCover_CONUS,
         relCoverA_NLTree = absNLTree_CONUS/totalAbsCover_CONUS,
         relCoverA_BLTree = absBLTree_CONUS/totalAbsCover_CONUS,
         relCoverA_C3Grass = absC3grass_CONUS/totalAbsCover_CONUS,
         relCoverA_C4Grass = absC4grass_CONUS/totalAbsCover_CONUS,
         relCoverA_Forb = absForb_CONUS/totalAbsCover_CONUS
         )

# make into a long format
predsOverTime_long <- predsOverTime %>% 
  select(year,x, y, spatialID, uniqueID, totalAbsCover_CONUS:relCoverA_Forb) %>% 
  pivot_longer(cols = c(totalAbsCover_CONUS:relCoverA_Forb)) %>% 
  mutate(year = as.integer(year))

# calculate the coefficient of variation for each cover type at each location
predsOverTime_cv <- predsOverTime_long %>% 
  group_by(uniqueID, name) %>% 
  summarize(cv = sd(value)/mean(value))
## add x,y 
predsOverTime_cv <- predsOverTime_cv %>% 
  left_join(predsOverTime %>% select(uniqueID, x, y))

# visualize change over time for a selection of points
predsOverTime_long %>% 
  filter(uniqueID %in% 1:50) %>% 
  filter(name %in% c( #"totalAbsCover_CONUS"#,  
    #"relCoverA_totalTree",   "relCoverA_totalHerb", 
    "relCoverA_shrub"    , "relCoverA_bareGround" , "relCoverA_NLTree",  "relCoverA_BLTree",       
"relCoverA_C3Grass" , "relCoverA_C4Grass",   "relCoverA_Forb" 
)) %>% 
  ggplot() + 
  facet_wrap(~uniqueID) + 
  geom_area(aes(x = year, y = value, fill = name)) + 
  # geom_line(data = predsOverTime_long %>% 
  #             filter(uniqueID %in% 1:50) %>% 
  #             filter(name %in% c("relCoverA_totalHerb", "relCoverA_totalTree")),
  #           aes(x = year, y = value, linetype = name), col = "black") +
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 90)) + 
  labs(title = "Change in modeled relative cover (all types sum to 1) by functional type from \n 2010 to 2020 for a subset of locations")

predsOverTime_cv %>%  
  ggplot() + 
  geom_boxplot(aes(x = name, y = cv, col = name)) +
  theme_minimal() + 
  labs(title = "Coefficients of variation of model-predicted  relative cover (all types sum to 1) \n at a given location for each functional group from 2010-2020",
        x = "functional group",
       y = "coefficient of variation") +
  theme(axis.text.x = element_text(angle = 90))

# show the locations of each point that is being used for this temporal analysis
#cropped_states <- cropped_states %>% 
 # rename(geometry = Shape)
 
ggplot() + 
  facet_wrap(~name) + 
  geom_sf(data = cropped_states, aes()) + 
  stat_summary2d(data = predsOverTime_cv, aes(x = x, y = y, z = cv), fun = ~ mean(.x, na.rm = TRUE)) +
  geom_hex(data = predsOverTime_cv, aes(x = x, y = y, col = cv), stat = "mean", na.rm = TRUE) + 
  guides(fill = guide_legend(title = "CV")) + 
  theme_minimal() + 
  labs(title = "Locations of points included in the random \n sample for examining temporal variation \n color-coded by the coefficient of variation")

Second option: relativize cover so that cover for all groups except trees sums to 1

## for contemporary data
contempPreds <- contempPreds %>% 
  mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS + 
                                           absBareGround_CONUS + absShrub_CONUS)) %>% 
  mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
         relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_NLTree = absNLTree_CONUS,
         relCoverB_BLTree = absBLTree_CONUS,
         relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS
         )

## for first climate model
predsFuture1 <- predsFuture1 %>% 
  mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS + 
                                           absBareGround_CONUS + absShrub_CONUS)) %>% 
  mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
         relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_NLTree = absNLTree_CONUS,
         relCoverB_BLTree = absBLTree_CONUS,
         relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS
         )

## for second climate model 
predsFuture2 <- predsFuture2 %>% 
  mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS + 
                                           absBareGround_CONUS + absShrub_CONUS)) %>% 
  mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
         relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_NLTree = absNLTree_CONUS,
         relCoverB_BLTree = absBLTree_CONUS,
         relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS
         )

Now, show the results in map form

# list of absolute cover names
responseNames <- c("relCoverB_totalTree", "relCoverB_totalHerb", "relCoverB_shrub", "relCoverB_bareGround",
                   "relCoverB_NLTree", "relCoverB_BLTree", "relCoverB_C3Grass", "relCoverB_C4Grass", "relCoverB_Forb",
                   "total_NonTree_AbsCover_CONUS")

relCover_B_Maps_contemp <- lapply(responseNames, 
                                      FUN = function(x) {
                                      
                                       ## contemporary absolute cover
                                        # rasterize response 
                                        resp_rast <- contempPreds %>% 
                                          terra::vect(geom = c("x", "y")) %>% 
                                          terra::set.crs(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE)   
                                        # get the extent of this particular raster, and crop it accordingly
                                        resp_rast_2 <- resp_rast %>% 
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        
                                        
                                         ## future model 1 absolute cover
                                        # rasterize response 
                                        resp_rast_future1 <- predsFuture1 %>% 
                                          terra::vect(geom = c("x", "y")) %>% 
                                          terra::set.crs(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE)   
                                        # get the extent of this particular raster, and crop it accordingly
                                        resp_rast_future1_2 <- resp_rast_future1 %>% 
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        # calculate differences between contemporary preds and future 1 preds (deltas)
                                        rast_deltas_future1 <- resp_rast_future1_2 - resp_rast_2
                                        
                                        ## future model 2 absolute cover
                                        # rasterize response 
                                        resp_rast_future2 <- predsFuture2 %>% 
                                          terra::vect(geom = c("x", "y")) %>% 
                                          terra::set.crs(crs(test_rast)) %>% 
                                          terra::rasterize(y = test_rast, 
                                                           field = x, 
                                                           fun = mean, na.rm = TRUE)   
                                        # get the extent of this particular raster, and crop it accordingly
                                        resp_rast_future2_2 <- resp_rast_future2 %>% 
                                          crop(ext(min(tempExt[,1]), max(tempExt[,1]),
                                                   min(tempExt[,2]), max(tempExt[,2])) 
                                               )
                                        # calculate differences between contemporary preds and future 1 preds (deltas)
                                        rast_deltas_future2 <- resp_rast_future2_2 - resp_rast_2
                                        
                                        
                                        if (x == "total_NonTree_AbsCover_CONUS") {
                                          # make a map of predictions w/ current data 
                                        resp_map <- ggplot() +
                                          geom_spatraster(data = resp_rast_2) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Predictions of Total Absolute Cover (non-tree) with contemporary data")) +
                                          scale_fill_gradient2(low = "lightblue",
                                                               #mid = "wheat" ,
                                                               high = "darkblue" , 
                                                               #midpoint = 0, 
                                                               limits = c(0,max(terra::values(resp_rast_2$mean), na.rm= TRUE)),  
                                                              na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist <- ggplot() + 
                                          geom_density(aes(contempPreds[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of Total Absolute Cover (non-tree) with contemporary data")) +
                                          theme_classic()
                                        
                                         # make a map of future preds w/ model 1
                                         
                                        resp_map_future1 <- ggplot() +
                                          geom_spatraster(data = rast_deltas_future1) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Future predictions of Total Absolute Cover (non-tree) \n compared to contemporary predictions (future - contemprary)  \n - climate model: BNU-ESM")) +
                                          scale_fill_gradient2(low = "orange",
                                                               mid = "white" ,
                                                               high = "purple" , 
                                                               midpoint = 0, limits = c(-2,2),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist_future1 <- ggplot() + 
                                          geom_density(aes(predsFuture1[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of relative cover of ",x, 
                                                      " \n in future scenario (model = BNU-ESM)")) +
                                          theme_classic()
                                        
                                        # make a map of future preds w/ model 2
                                        resp_map_future2 <- ggplot() +
                                          geom_spatraster(data = rast_deltas_future2) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Future predictions of Cover (non-tree), \n compared to contemporary predictions (future - contemprary) \n- climate model: IPSL-CM5A-MR (France)")) +
                                          scale_fill_gradient2(low = "orange",
                                                               mid = "white" ,
                                                               high = "purple" , 
                                                               midpoint = 0, limits = c(-2,2),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist_future2 <- ggplot() + 
                                          geom_density(aes(predsFuture2[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of relative cover of ",x, 
                                                      " \n in future scenario (model = IPSL-CM5A-MR (France))")) +
                                          theme_classic()
                                        
                                        } else {
                                          # make a map of predictions w/ current data 
                                        resp_map <- ggplot() +
                                          geom_spatraster(data = resp_rast_2) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Predictions of relative cover of ",x,
                                                                        ", \n when all functional groups except trees sum to 1")) +
                                          scale_fill_gradient2(low = "brown",
                                                               mid = "wheat" ,
                                                               high = "darkgreen" , 
                                                               midpoint = 0, limits = c(0,1),  
                                                              na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist <- ggplot() + 
                                          geom_density(aes(contempPreds[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of relative cover of ",x)) +
                                          theme_classic()
                                        
                                         # make a map of future preds w/ model 1
                                         
                                        resp_map_future1 <- ggplot() +
                                          geom_spatraster(data = rast_deltas_future1) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Future predictions of relative cover of ",x,", \n compared to contemporary predictions (future - contemprary)  \n - climate model: BNU-ESM")) +
                                          scale_fill_gradient2(low = "orange",
                                                               mid = "white" ,
                                                               high = "purple" , 
                                                               midpoint = 0, limits = c(-1,1),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist_future1 <- ggplot() + 
                                          geom_density(aes(predsFuture1[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of relative cover of ",x, 
                                                      " \n in future scenario (model = BNU-ESM)")) +
                                          theme_classic()
                                        
                                        # make a map of future preds w/ model 2
                                        resp_map_future2 <- ggplot() +
                                          geom_spatraster(data = rast_deltas_future2) +
                                          geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
                                                    st_crop(st_bbox(plotObs_2)),fill=NA ) +
                                                    labs(title = paste0("Future predictions of relative cover of ",
                                                                        x,
                                                                        ", \n compared to contemporary predictions (future - contemprary) \n- climate model: IPSL-CM5A-MR (France)")) +
                                          scale_fill_gradient2(low = "orange",
                                                               mid = "white" ,
                                                               high = "purple" , 
                                                               midpoint = 0, limits = c(-1,1),  na.value = "lightgrey") + 
                                          xlim(st_bbox(plotObs_2)[c(1,3)]) + 
                                          ylim(st_bbox(plotObs_2)[c(2,4)]) + 
                                          theme_classic()
                                        
                                        # make a histogram
                                        resp_hist_future2 <- ggplot() + 
                                          geom_density(aes(predsFuture2[,x]), col = "black", fill = "darkgrey") + 
                                          xlab(paste0("Predictions of relative cover of ",x, 
                                                      " \n in future scenario (model = IPSL-CM5A-MR (France))")) +
                                          theme_classic()
                                        }
                                        
                                        
                                        return(list("rast" = resp_rast_2, 
                                                    "map" = resp_map, 
                                                    "hist" = resp_hist,
                                                    "rast_future1" = resp_rast_future1_2,
                                                    "map_future1" = resp_map_future1, 
                                                    "hist_future1" = resp_hist_future1,
                                                    "rast_future2" = resp_rast_future2_2,
                                                    "map_future2" = resp_map_future2, 
                                                    "hist_future2" = resp_hist_future2))
                                        })


names(relCover_B_Maps_contemp) <- c("relCoverB_totalTree", "relCoverB_totalHerb", "relCoverB_shrub", "relCoverB_bareGround",
                   "relCoverB_NLTree", "relCoverB_BLTree", "relCoverB_C3Grass", "relCoverB_C4Grass", "relCoverB_Forb",
                   "total_NonTree_AbsCover_CONUS")
ggarrange(
ggarrange(relCover_B_Maps_contemp$relCoverB_totalTree$map, 
          relCover_B_Maps_contemp$relCoverB_totalTree$map_future1,
          relCover_B_Maps_contemp$relCoverB_totalTree$map_future2,
          relCover_B_Maps_contemp$relCoverB_totalTree$hist,
           relCover_B_Maps_contemp$relCoverB_totalTree$hist_future1, 
          relCover_B_Maps_contemp$relCoverB_totalTree$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_B_Maps_contemp$relCoverB_totalHerb$map, 
          relCover_B_Maps_contemp$relCoverB_totalHerb$map_future1,
          relCover_B_Maps_contemp$relCoverB_totalHerb$map_future2,
          relCover_B_Maps_contemp$relCoverB_totalHerb$hist,
          relCover_B_Maps_contemp$relCoverB_totalHerb$hist_future1, 
          relCover_B_Maps_contemp$relCoverB_totalHerb$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_B_Maps_contemp$relCoverB_shrub$map, 
          relCover_B_Maps_contemp$relCoverB_shrub$map_future1,
          relCover_B_Maps_contemp$relCoverB_shrub$map_future2,
          relCover_B_Maps_contemp$relCoverB_shrub$hist,
          relCover_B_Maps_contemp$relCoverB_shrub$hist_future1, 
          relCover_B_Maps_contemp$relCoverB_shrub$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_B_Maps_contemp$relCoverB_bareGround$map, 
          relCover_B_Maps_contemp$relCoverB_bareGround$map_future1,
          relCover_B_Maps_contemp$relCoverB_bareGround$map_future2,
          relCover_B_Maps_contemp$relCoverB_bareGround$hist,
          relCover_B_Maps_contemp$relCoverB_bareGround$hist_future1, 
          relCover_B_Maps_contemp$relCoverB_bareGround$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_B_Maps_contemp$relCoverB_NLTree$map, 
          relCover_B_Maps_contemp$relCoverB_NLTree$map_future1,
          relCover_B_Maps_contemp$relCoverB_NLTree$map_future2,
          relCover_B_Maps_contemp$relCoverB_NLTree$hist,
          relCover_B_Maps_contemp$relCoverB_NLTree$hist_future1, 
          relCover_B_Maps_contemp$relCoverB_NLTree$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_B_Maps_contemp$relCoverB_BLTree$map, 
          relCover_B_Maps_contemp$relCoverB_BLTree$map_future1,
          relCover_B_Maps_contemp$relCoverB_BLTree$map_future2,
          relCover_B_Maps_contemp$relCoverB_BLTree$hist,
          relCover_B_Maps_contemp$relCoverB_BLTree$hist_future1, 
          relCover_B_Maps_contemp$relCoverB_BLTree$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_B_Maps_contemp$relCoverB_C3Grass$map, 
          relCover_B_Maps_contemp$relCoverB_C3Grass$map_future1,
          relCover_B_Maps_contemp$relCoverB_C3Grass$map_future2,
          relCover_B_Maps_contemp$relCoverB_C3Grass$hist,
          relCover_B_Maps_contemp$relCoverB_C3Grass$hist_future1, 
          relCover_B_Maps_contemp$relCoverB_C3Grass$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_B_Maps_contemp$relCoverB_C4Grass$map, 
          relCover_B_Maps_contemp$relCoverB_C4Grass$map_future1,
          relCover_B_Maps_contemp$relCoverB_C4Grass$map_future2,
          relCover_B_Maps_contemp$relCoverB_C4Grass$hist,
          relCover_B_Maps_contemp$relCoverB_C4Grass$hist_future1, 
          relCover_B_Maps_contemp$relCoverB_C4Grass$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ggarrange(relCover_B_Maps_contemp$relCoverB_Forb$map, 
          relCover_B_Maps_contemp$relCoverB_Forb$map_future1,
          relCover_B_Maps_contemp$relCoverB_Forb$map_future2,
          relCover_B_Maps_contemp$relCoverB_Forb$hist,
          relCover_B_Maps_contemp$relCoverB_Forb$hist_future1, 
          relCover_B_Maps_contemp$relCoverB_Forb$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ), 
ggarrange(relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$map, 
          relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$map_future1,
          relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$map_future2,
          relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$hist,
          relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$hist_future1, 
          relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$hist_future2,
          ncol = 3,
          nrow = 2,
          heights = c(1,.35)
          ),
ncol = 1
)

Now, show how the relativized cover values change over time

# calculate relativized cover (version w/ all functional groups summing to 1)
predsOverTime <- predsOverTime %>% 
mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS + 
                                         absBareGround_CONUS + absShrub_CONUS)) %>% 
  mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
         relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_NLTree = absNLTree_CONUS,
         relCoverB_BLTree = absBLTree_CONUS,
         relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
         relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS
         )

# make into a long format
predsOverTime_long <- predsOverTime %>% 
  select(year,x, y, spatialID, uniqueID, total_NonTree_AbsCover_CONUS:relCoverB_Forb) %>% 
  pivot_longer(cols = c(total_NonTree_AbsCover_CONUS:relCoverB_Forb)) %>% 
  mutate(year = as.integer(year))

# calculate the coefficient of variation for each cover type at each location
predsOverTime_cv <- predsOverTime_long %>% 
  group_by(uniqueID, name) %>% 
  summarize(cv = sd(value)/mean(value))
## add x,y 
predsOverTime_cv <- predsOverTime_cv %>% 
  left_join(predsOverTime %>% select(uniqueID, x, y))

# visualize change over time for a selection of points
predsOverTime_long %>% 
  filter(uniqueID %in% 1:50) %>% 
  filter(name %in% c( #"total_NonTree_AbsCover_CONUS",  "relCoverB_totalTree",    
#"relCoverB_totalHerb", 
    "relCoverB_shrub"    , "relCoverB_bareGround" , "relCoverB_NLTree",  "relCoverB_BLTree",       
"relCoverB_C3Grass" , "relCoverB_C4Grass",   "relCoverB_Forb" )) %>% 
  ggplot() + 
  facet_wrap(~uniqueID) + 
  geom_area(aes(x = year, y = value, fill = name)) + 
  geom_line(data = predsOverTime_long %>% 
              filter(uniqueID %in% 1:50) %>% 
              filter(name %in% c("relCoverB_totalHerb", "relCoverB_totalTree")),
            aes(x = year, y = value, linetype = name), col = "black") +
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 90)) + 
  labs(title = "Change in modeled relative cover (all Non-tree types sum to 1) by functional type from \n 2010 to 2020 for a subset of locations")

predsOverTime_cv %>%  
  ggplot() + 
  geom_boxplot(aes(x = name, y = cv, col = name)) +
  theme_minimal() + 
  labs(title = "Coefficients of variation of model-predicted  relative cover (all Non-tree types sum to 1) \n at a given location for each functional group from 2010-2020",
        x = "functional group",
       y = "coefficient of variation") +
  theme(axis.text.x = element_text(angle = 90))

# show the locations of each point that is being used for this temporal analysis
#cropped_states <- cropped_states %>% 
 # rename(geometry = Shape)
 
ggplot() + 
  facet_wrap(~name) + 
  geom_sf(data = cropped_states, aes()) + 
  stat_summary2d(data = predsOverTime_cv, aes(x = x, y = y, z = cv), fun = ~ mean(.x, na.rm = TRUE)) +
  geom_hex(data = predsOverTime_cv, aes(x = x, y = y, col = cv), stat = "mean", na.rm = TRUE) + 
  guides(fill = guide_legend(title = "CV")) + 
  theme_minimal() + 
  labs(title = "Locations of points included in the random \n sample for examining temporal variation \n color-coded by the coefficient of variation")